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ABSTRAGT 

We revisit the problem of low-mass pre-main-sequence (PMS) stellar evolution and its observational 
consequences for where stars fall on the Hertzsprung- Russell diagram (HRD). In contrast to most 
previous work, our models follow stars as they grow from small masses via accretion, and we perform 
a systematic study of how the stars' HRD evolution is influenced by their initial radius, by the radiative 
properties of the accretion flow, and by the accretion history, using both simple idealized accretion 
histories and histories taken from numerical simulations of star cluster formation. We compare our 
numerical results to both non- accreting isochrones and to the positions of observed stars in the HRD, 
with a goal of determining whether both the absolute ages and the age dispersions inferred from non- 
accreting isochrones are reliable. We show that non-accreting isochrones can sometimes overestimate 
stellar ages for more massive stars (those with effective temperatures above ^ 3500 K), thereby 
explaining why non-accreting isochrones often suggest a systematic age difference between more and 
less massive stars in the same cluster. However, we also find the only way to produce a similar 
overestimate for the ages of cooler stars is if these stars grow from ~ 0.01 Mq seed protostars that are 
an order of magnitude smaller than predicted by current theoretical models, and if the size of the seed 
protostar correlates systematically with the final stellar mass at the end of accretion. We therefore 
conclude that, unless both of these conditions are met, inferred ages and age spreads for cool stars 
are reliable, at least to the extent that the observed bolometric luminosities and temperatures are 
accurate. Finally, we note that the time-dependence of the mass accretion rate has remarkably little 
effect on low-mass stars' evolution on the HRD, and that such time-dependence may be neglected for 
all stars except those with effective temperatures above ~ 4000 K. 

Subject headings: accretion — Hertzsprung-Russell diagram — stars: evolution — stars: formation 
— stars: pre-main sequence 



1. INTRODUCTION 

Pre-main-sequence (PMS) stars in low-mass star form- 
ing regions show a sizable luminosity spread when 
place d on the Hertzsp rung-Russell diagram (HRD) 
(e.g., iHillenbrandl 120091) . This spread translates into 
a significant dispersion in inferred stellar ages that 
recor ds the past star formation ac t ivity in each region 
(e.g., i D'Antona fc MazzitcUi' Togi; 'Baraffc ct al. 1998': 
Siess et al. 2000: PaUa & Stabler 2000, 1999; Hartinann 
2001L l2003h . However, the idea that star clusters form 
over an extended period is subject to extensive de- 
bate on both observational and theoretical grounds (e.g. , 
ElmcR-rccn 2000: Hartmann et al."2001': ' Tan et al.l l2006l: 
Krumholz & Tan 2007; Evans et al. 2009^, and several 
authors have claimed that the dispersion of stellar lu- 
minosities does not reflect a real age spread. Mem- 
bers of young binaries and multiples exhibit a tighter 
age correlation, supporting the existence of an intrin- 
sic age distribution. However, a luminosity spread per- 
sists even among such systems, and some companions 
display a substantial age mismatch (jPrato et al.l 120031 : 
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iStassun et al|[2008t iKraus fc Hillenbrandll2009[ ). Deriv- 
ing stellar ages is complicated from an observational 
standpoint. Galculation of stellar bolometric luminosi- 
ties is beset by uncertainties in extinction, photomet- 
ric variability, and unresolved multiplicity. Galibration 
between the stellar spectral type and effective tempera- 
ture is also not trivial. In some cases, it can be demon- 
strated that observational uncertainties alone are suffi- 
cient to induce an age spre ad of > 10 Myr an d mask 
a coeval stellar population (' Slesnick et al.l [20081 ). How- 
ever, [Da Rio et al. (2010a. b) carefully model these un- 
certainties and conclude that these effects alone cannot 
reproduce the entire spre ad. Other age indica tors such 
as stella r rotation rate (iLittlefair et all I2O1O0 . surface 
gravitv (iSle snick et al.' '2008"), and lithium abundances 
(jSestito et al. 2008) also support the idea that the in- 
ferred age spreads are real, but are each subject to sig- 
nificant challenges. 

Apart from the observational uncertainties, physical 
mechanisms may be responsible for a portion of the ob- 
served HRD scatter. For the purpose of inferring stel- 
lar ages, it is usually assumed that PMS stars first ap- 
pear along a "bi rthline" in the HRD when mass accre- 
tion ceases (e.g., IPalla fc Stahleii [19901 : IHartmann et al.l 
Il997f ). However, luminosities of younger embedded stars 
(Glass and I sources) that are presumably still accreting 
also show a wide spread, and a fraction of them have lu- 
minosities much lower than the values ex pected from the 
standard birthline fe.g.. iKenvon et al.lll990tlEvans et al.l 
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l2009bHEnoch et al.|[2009l ). A solution for this "luminos- 
ity problem" is the scenario that mass accretion takes 
place very time-dependently, repeating burst-like accre- 
tion phases and quiescent phases. Recent numerical 
simulations suggest that such episodic mass accretion is 
caused by gravitational fragmentation of a circmustella r 
disk fe-g- lVorobvov fc Basull2005HMachida et al.ll2QTTI ). 
thou gh radiative warm ing from protostars alleviates it 
(e.g.. lOffner et al.l l2009:l. Regardless of the ultimate ex- 
planation for the luminosities of Class and I sources, 
the existence of young stars that fall well below the puta- 
tive birth line is strong evidence that we must extend our 
P MS evolution model s to include the accretion phase. 

iBaraffe et all (|2009l hereafter BCG09) study protostel- 
lar evolution with various episodic mass accretion histo- 
ries and examine the resultant spread of PMS stars in 
the HRD. They argue that PMS stars of the same mass 
and age show some scatter in the HRD owing to varia- 
tion of the early evolution resulting from complex accre- 
tion histories. However, BCG09 simultaneously vary not 
only their accretion histories, but also their initial stel- 
lar models and the radiative properties of the accretion 
flow. Because they change these parameters in correlated 
ways and do not perform a systematic survey of parame- 
ter space, it is not clear which of these effects drives their 
results. Nor is it clear whether the results they generate 
via their parameter choices are consistent with observed 
HRDs of clusters. Consequently, it is still unclear how 
much vigorous time-dependent accretion histories influ- 
ence protostellar evolution. 

In this paper, we aim to resolve this question by per- 
forming a systematic study of how PMS evolutionary 
tracks change as we alter the accretion history, the ini- 
tial models, and the thermal efficiencies of mass accre- 
tion. We perform a systematic survey of parameter space 
in order to understand how each of these factors affects 
protostellar evolution. This enables us to answer the 
question of whether variation in any of these quantities 
could produce the appearance of an age spread in a pop- 
ulation that is actually coeval. 

The structure of the paper is as follows. In Section 
[21 we briefly explain our numerical method for modeling 
protostellar evolution. Section [3] is the main part of the 
paper, where the numerical results are presented. First, 
we investigate how different accretion histories influence 
protostellar evolution in Section 13.11 We next investi- 
gate protostellar evolution with differing initial models 
in Section 13.21 and with differing thermal efficiencies in 
13.31 In Section 2] we combine all these results to draw 
general conclusions about the reliability of age and age 
spread estimates from PMS evolutionary tracks. Section 
[5] contains the summary and discussion. 

2. PROTOSTELLAR EVOLUTION MODELS 

We model protostel lar evolution using the nu meri- 
cal code describ e d by iHosokawa fc Omukail (|2009[ ) and 
iHosokawa et al.l ()2010D . The code numerically solves 
the four stellar structure equations, taking into ac- 
count mass accretion. For the follo wing calculations, 
we a dopt the OPAL opacity tables (jlglesias fc Rogera 
I1996D for high tempera ture T > 7000 K, and o t her ta bles 
based on the work by [Alexander fc FergusonI ()1994[ ) for 
the lower temperature. We employ mixing-length the- 
ory for heat transport in convective layers with a con- 



stant ratio of the mixing length to the pressure scale 
height of 1.5. We confirmed that our code reproduces 
the ca lculati ons bylS tahler et al.' Cl980D. |PaIla fc Stahleil 
(I1990D . and iPalla fc Stabler (1992) in both the lim- 
its of hot spherical accretion and cold disk accretion, 
w hich we explain in mo re detail bel ow (see append ixes 
in IHosokawa fc Omukai 2009: Hosok awa et aLlboiOO . 

We refer the reader to the Hosokawa et al. papers for 
full details of the numerical method, but one parameter 
is particularly important for the results of this paper. 
The thermal efficiency of mass accretion, i.e., the entropy 
carried into the star with accreting material, is a key 
parameter for protostellar evolution. Since a protostar 
grows by accretion, the average entropy in the stellar 
interior becomes higher with higher thermal efficiency. 
For a star with a fixed mass, the stellar radius is larger 
for higher interior entropy content. Thus, we naively 
expect that, even for fixed accretion history, protostars 
will have larger radii if the accretion flows onto them 
have higher thermal efficiencies. 

Despite efforts in previous work, however, the concrete 
value of the thermal efficiency in low-mass star forma- 
tion is not well-constrained. Here we address this un- 
certainty by considering two limiting cases, representing 
cartoon versions of two different accretion fiow geome- 
tries: "hot" spherical accretion, and "cold" disk accre- 
tion. In the hot accretion case, we envision that an ac- 
cretion flow directly hits the stellar surface and forms 
an accretion shock front. The accretion flow may ar- 
rive in a disk, but in the hot case we imagine that the 
disk is thick enough so that the accretion column covers 
much of the stellar surface. As a result, a small frac- 
tion of the heat generated at the shock front is carried 
into the stellar interior. In this limit we solve for the 
steady structure of the gas accretion envelope as well as 
the stellar interior; the two are connected at the stel- 
lar surface with accretion shock jump conditi ons (e.g., 
iStahler et all 119801: IHosokawa fc Omukail 120091 ). In con- 
trast, in the cold accretion case we envision that accret- 
ing gas initially falls onto a circumstellar disk and then 
reaches the stellar surface via a thin accretion column 
connecting the disk and star. As a result, most of the 
stellar photosphere is not covered by an accretion column 
and is able to radiate freely. Accreting gas softly settles 
on the stellar surface, and when it is incorporated into 
the star it has the same entropy as gas in the stellar pho- 
tosphere. In this case we do not solve for the structure 
of the accretion flow, and we instead ad opt the ordinary 
photospheric boundary con dition (e.g., IPalla fc Stahleii 
fl99l IHosokawa et al.ll20Tol) . 

Our treatment of boundary conditions differs slightly 
from that of BCG09, who modeled the thermal efficiency 
with a parameter a, the fraction of accreting internal en- 
ergy absorbed by the star. However, our limiting cases 
of hot and cold accretion just correspond to their a = 1 
and cases respectively. The only other difference be- 
tween our and BCG09's method is that BCG09 assume 
instantaneous and uniform mixing of accreting material 
in the stellar interior (Sicss & Forcstini 1996). In this 
case, the entropy of newly accreted material is assumed 
to be the same as the local values in the stellar interior, 
whereas in our cold case it is assumed to match the stel- 
lar photosphere. Note that, even in the hot case, only 
a small fraction of the accretion energy goes into heat- 
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Fig. 1. — Time-dependent accretion histories adopted for cal- 
culations of protostellar evolution. The thick solid line presents 
the sample accretion history tak en from the numeri cal simulations 
of low-mass star formation by ' Offner et al.l 1 120091) . and used in 
model mO-C. The thin solid line shows the vigorous episodic accre- 
tion case used in model mE-C, where a burst-like accretion phase 
at 10~* Mq yr~^ over 100 years is interspersed with a quiescent 
phase at 3 X 10"'' Mq yr~^ over 1000 years. The accretion history 
over initial 40000 years is shown for this case. The small window 
enlarges the evolution over the initial 4000 years. 



ing the stellar matter. This is similar to the a — 1 case 
in BCG09, where the accretion energy acts as a uniform 
heating source that is dist ributed uniforml y throughout 
the stellar interior (e.g., see lSiess et al.|[l997[ ). As a result 
of this treatment, most of this energy escapes from the 
star without being absorbed by the stellar matter. Thus, 
the term "hot mass accretion," in both our treatment and 
BCG09's, indicates that mass accretion increases the av- 
erage entropy in the stellar interior, not that the accre- 
tion flow is completely radiatively inefficient. 

3. RESULTS 

3.1. Accretion History Variation 

First, we examine how variations in protostellar ac- 
cretion histories affect stars' HRD evolution. To this 
end, we calculate a series of models, summarized in Ta- 
ble [1^. In order to isolate dependence on the accretion 
history from other effects, all these calculations use the 
same initial model and boundary conditions. This ini- 
tial model co nsists of a 0.0 1 Mp, star with radiative 
interior, as in iStahler et al.l (|1980l ). and an initial ra- 
dius of 1.5 Rq. Note that the value of 1.5 Rq is a bit 
sm aller than the radius o f the s eed protostar calculated 
bv lMasunaga fc Inutsukal (|200a hereafter MIOO), 4 Rq. 
The boundary condition for all these models is cold 
accretion. We note that BCG09 obtained protostellar 
evolutionary tracks that deviate from the non-accreting 
isochrones substantially only for their cold cases (their 
a = 0), and this motivates us to focus on cold accretion 
first. 

3.1.1. Episodic Accretion 

We first consider five distinct accretion histories with 
varying degrees of episodic variation, but with only fac- 
tor of 2 level changes in the average accretion rate. 
This enables us to explore the sensitivity of stars' HRD 
evolution to the level of variability in their accretion his- 
tories. Model mO-C uses an accretion history taken 
fro m n umerical simulations of low-mass star formation 
by lOffner et al.. (,2009i) and is illustrated in Figure [TJ 
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Fig. 2. — Evolution of the stellar radius versus stellar mass for 
cases mC5-C (red), mE-C (magenta), mOx0.5-C (blue solid), mO- 
C (blue dashed), and mOx2-C (blue dotted), among which the 
accretion histories differ but the initial and boundary conditions 
do not. 



In this case, the accretion rate gradually decreases over 
~ 0.1 Myr, and the stellar mass finally reaches 0.45 Mq. 
Models mOx0.5-C and mOx2-C use the same accretion 
history, but they are scaled by factors of 0.5 and 2, re- 
spectively, to give final masses of 0.23 Mq and 0.9 Mq. 
In contrast, model mC5-C uses a fixed accretion rate of 
M = 10"^ Mq yr-^ We use this model to produce 0.23, 
0.45, and 0.9 Mq stars, as in the niO-C models, simply 
by stopping accretion once the star has reached the de- 
sired mass. Finally, model mE-C (also shown in Figure 
[T]) represents an extreme case of variability: an episodic 
mass accretion history where burst-like accretion events 
at M = lO"** Mq yr~^ over 100 years are interspersed 
with quiescent phases of accretion at 3 x 10~^ Mq yr~^ 
over 1000 years. This model is similar to the episodic 
accretion histories formulated by BCG09. As with mC5- 
C, we use this model to produce 0.23 Mq, 0.45 Mq, and 
0.9 Mq stars simply by turning off mass accretion once 
the stellar mass reaches the target value. 

Figure[2]presents the evolution of the stellar radius un- 
til mass accretion ceases in each case. We see that the ba- 
sic evolution is similar for all cases. The stellar radius ini- 
tially decreases with increasing stellar mass. The temper- 
ature in the stellar interior rises during this initial con- 
traction. When the stellar mass reaches M^ ~ 0.07 Mq, 
deuterium burning begins and the stellar interior be- 
comes fully convective. For some time after the ignition 
of deuterium, temperature at the stellar center remains 
constant at ~ 10® K due to the very strong temperature- 
dependence of the deuterium burning rate. This is the 
so-called ther mostat effect of deuterium burning (e.g., 
i StahledfTgsl . The stehar radius mcreases m propor- 
tion to the stellar mass during this phase. The deu- 
terium concentration in the stellar interior significantly 
decreases with increasing the stellar mass dStahlcr 198| 
iHartmann et al.lll997[ ). Finally at > 0.2 Mq, the 
deuterium concentration is so low that the thermostat 
effect becomes inoperative and the central temperature 
increases again. Variation of the accretion histories only 
slightly influences the evolution of the stellar radius. 

Figure shows the stellar positions in the HRD at 

^ Note that this figure is intended to facilitate comparison be- 
tween the models and data for PMS stars, for which mass accretion 
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TABLE 1 

Model parameters and results 



Case Accretion History Boundary Condition -R*,o (^©) A<f*,d (Mq) (Mq) R*j (Rq) (kyr) 



la. Fixed Initial and Boundary Conditions, Varying Accretion History 



mC5-C 


10-5 Mq yr-i 


C 


1.5 


0.074 


0.9 


1.3 


90 


niE-C 


episodic 


C 


1.5 


0.07 


0.9 


1.8 


90 


mO-C 


simulation'' 


c 


1.5 


0.075 


0.45 


1.3 


110 


mOx2-C 


simulation'' 


c 


1.5 


0.077 


0.9 


1.1 


110 


mOx0.5-C 


simulation'' 


c 


1.5 


0.073 


0.23 


1.3 


110 


mC4-C 


10-4 Mq yr-1 


c 


1.5 


0.076 


0.9 


0.92 


9 


mC6-C 


10-6 Mq yr-i 


c 


1.5 


0.1 


0.9 


1.4 


900 




lb. Fixed Accretion 


History and Boundary Conditions, Varying Initial Conditions 






mC5-C-Ri-8 


10-5 Mq yr-1 


C 


8.0 


0.093 


0.9 


1.8 


90 


mC5-C-Ri3.7 


10-5 Mq yr-i 


C 


3.7 


0.09 


0.9 


1.7 


90 


mC5-C 


10-5 Mq yr-i 


C 


1.5 


0.074 


0.9 


1.3 


90 


mC5-C-Ril 


10-5 Mq yr-i 


c 


1.0 


0.06 


0.9 


1.1 


90 


mC5-C-Ri0.65 


10-5 Mq yr-i 


c 


0.65 


0.049 


0.9 


0.85 


90 


mC5-C-Ri0.3 


10-5 Mq yr-i 


c 


0.3 


0.035 


0.9 


0.64 


90 


mC5-C-Ri0.25 


10-5 Mq yr-l 


c 


0.25 


0.033 


0.9 


0.4 


90 


mC5-C-Ri0.2 


10-5 Mq yr-1 


c 


0.2 


0.033 


0.9 


0.29 


90 



Ic. Fixed Accretion History and Initial Conditions, Varying Boundary Conditions 



mC5-C-Ri3.7 


10" 


'5 


Mq yr- 


-1 


C 


3.7 


0.09 


0.9 


1.7 


90 


mC5-C 


10- 


-5 


Mq yr- 


-1 


C 


1.5 


0.074 


0.9 


1.3 


90 


mC5-H 


10- 


-5 


Mq yr- 


-1 


H 


3.7 


0.34 


0.9 


4.6 


90 


mC5-HC0.3 


10- 


-5 


Mq yr- 


■ 1 


H -i> C (0.3 Mq^ 


3.7 (3.0)'' 


0.33 


0.9 


4.3 


90 


mC5-HC0.1 


10- 


-5 


Mq yr- 


■ 1 


H C (0.1 Mq)'' 


3.7 (2.6)'' 


0.25 


0.9 


3.6 


90 


mC5-HC0.03 


10- 


-5 


Mq yr- 


-1 


H ^ C (0.03 Mq)^ 


3.7 (3.2)'' 


0.17 


0.9 


2.9 


90 



Note. — Col. 2: see main text for details of the accretion histories wo use; Col. 3: H — hot accretion, C — cold accretion; Col. 
4: initial stellar radius when M*: — 0.01 Mq; Col. 5: stellar mass when core deuterium burning begins; Col. 6: final stellar mass for 
the most massive star we produce with these conditions; for most cases we also run to a series of smaller final masses; Col. 7: stellar 
radius at the end of accretion for the most massive case wo run; Col. 8: time when accretion ends for the most massive case we run 
^ For run mO-C the accretion history is taken from the simulations of lOffner et al.] | |2009[) (see text for details). Runs mOx0.5-C and 
mOx2-C use the same accretion history, scaled by factors of 0.5 and 2, respectively. 

^ For run mC5-HCa:, the boundary condition is switched from hot to cold once the stellar mass reaches x Mq. The quantities given 
in parentheses are the stellar mass and radius when this switch occurs. 



t ^ 0.3 Myr, 1 Myr, 3 Myr, and 10 Myr after mass 
accretion begins for each mass accretion history. The 
snapshot at t — 0.3 Myr shows the steUar positions just 
after mass accretion ceases. Reflecting the minimal vari- 
ation in the stars' radii shown in Figure [21 the stellar 
positions only show a small spread. This means that 
the concept of the birthline is valid even with variable 
accretion histories in the limiting case of cold mass ac- 
cretion, provided the initial state is the same from star 
to star. At t > 0.3 Myr, the stars gradu ally approach 
the ZAMS line, descending in the HRD ([Henvev et alj 
[19551: IHavashi|[l96ll : IHavashi fc Nakandll963D . The snap- 
shot at t = 1 Myr clearly shows that the stars are below 
the 1 Myr isochrone for non-accreting protostars. This 
offset is larger for the higher-mass stars. In particular, 
the 0.9 Mq stars are close to the isochrone of 10 Myr. 
The offset decreases with time, but still remains even at 
t = 10 Myr for 0.9 Mq stars. 

This divergence between the accreting evolutionary 
histories and the non-accreting isochrones is easy to un- 
derstand. The isochrones for non-accreting protostars 
are derived assuming a large initial radius of i?*,o ~ 

has presumably ceased. In this and all subsequent figures we omit 
the accretion luminosity for the early evolutionary tracks. If the 
accretion luminosity is added, the tracks shift upward in the HRD 
and no longer pass through the data. 



10 Rq. The model PMS stars then contract from this ini- 
tial state by radiating away their energy, reaching smaller 
radii at larger stellar ages. On the other hand, with ther- 
mally inefficient accretion the stellar radius remains as 
small as i?* < 1 Rq during mass accretion. As a re- 
sult, stellar ages are overestimated using the isochrones 
for non-accreting protostars in such cases. However, we 
stress that this effect has nothing to do with the time 
variability of the mass accretion rate. 

3.1.2. Varying Mean Accretion Rates 

Next we examine protostellar evolution over a greater 
range of mass accretion rates but with no time variabil- 
ity. In addition to case mC5-C, we consider cases mC4- 
C and mC6-C, which have constant accretion rates of 
10~* Mq yr~^ and 10~^ AIq yr~^, respectively. Figure 
[31 is the same as Figure [31 but for these cases. We see 
larger variation of the tracks among these cases than in 
Figure [31 but the level of variation is still significantly 
smaller than the observed range of data, particularly for 
the lowest mass and effective temperature. There is little 
difference between cases mC5-C and mC6-C, where ac- 
cretion rates take typical values for low-mass star forma- 
tion of M < 10^5 Mq yr-i. The mC4-C track slightly 
deviates from these tracks. In this case, however, the 
0.9 Mq star always lies below the 10 Myr isochrone, far 
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Fig. 3. — Snapshots of stellar positions in the HRD for varying accretion histories with the same initial and boundary conditions. The 
four panels show the snapshots at times 0.3 Myr (upper left), 1 Myr (upper right), 3 Myr (lower left), and 10 Myr (lower right) after the 
start of mass accretion. The evolutionary tracks until that time are also plotted in the panels. (Note that the luminosity plotted here 
is only the stellar luminosity; accretion luminosity is not included.) The different colors represent different accretion histories: constant 
accretion at 10~^ Mq yr~^ (red, case mC5-C), episodic accretion (magenta, case mE-C), and decreasing accretion (blue, cases mO-C, 
mOx2-C, mOx0.5-C). The input parameters in each case are summarized in Table [l^. In each panel, the symbols mark the positions of 
stars whose masses are 0.9 M p, (asterisks), 0.45 Mq (pluses), and 0.23 Mq (crosses). The thick dot-solid line indicates the positions of 
ZAM S stars fSiess et al.' 2000). The dashed lines represent the isochrones o f 1 Myr and 10 My r for non-accreting protostars ( Baraf fe et aLl 
I199SI) . Observational data is taken from Gatti et al. (2006, open squares). IGatti et al.l Il2008i . filled squares). ;Muzerolle et al. (200 51. open 
triangles), and iPeterson et al.i ([200^, filled triangles). 



from the locations where observed protostars lie. We 
therefore conclude that purely cold accretion at rates of 
10~* Mq yr~^, while not physically forbidden, docs not 
appear to actually occur in observed star clusters, at least 
for 0.9 Mq stars. The problem that many evolution- 
ary scenarios involving cold accretion overpopulate the 
region at low L and high Toff is one we will encounter 
repeatedly in the rest of this paper. 

In summary, we conclude that star-to-star differences 
in either the overall accretion rate or the degree to which 
the accretion rate varies in time have a very limited effect 
on protostellar evolution. As illustrated by Figures [3] 
and m varying accretion histories can cause absolute age 
estimates to be wrong for more massive stars, but cannot 



explain the observed broad spread of PMS stars in the 
HRD. If PMS stars' initial state and accretion boundary 
conditions were fixed, then a co-eval population would 
form a much tighter sequence in the HRD than what we 
actually observe, even if their accretion histories varied 
wildly. 

3.2. Initial Model Variation 

We next study protostellar evolution by varying the 
initial model, while fixing the boundary condition to 
cold mass accretion and the accretion rate to M = 
10^^ Mq yr~^. We choose cold accretion here, be- 
cause protostellar evolution is only sensitive to the initial 
model in the cold case, not the hot one (,Stahler.,1988l: 
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Fig. 4. — Same as Figure \3\ but for different constant accretion rates. The different colors donate differences of the accretion rates: 
10~^ Mq yr~^ (green, case mC6-C), 10~^ Mq yr~^ (red, case mC5-C), and 10~* Mq yr~^ (blue, case mC4-C). Parameters in each case 
are summarized in Table [Tk. 



iHartmann et al.lll997l ) - we defer a detailed discussion of 
this issue to Section [331 Throughout these calculations, 
we assume a fixed initial deuteri um abu n dance , which in 
principle may vary somewhat. iStahleii ()1988D explored 
the effect of varying deuterium abundance on the stellar 
radius. He showed that larger amounts of deuterium led 
to somewhat larger stars. However, this had only a small 
effect on the total thermal energy due to the thermostat 
effect of deuterium burning, which acts to regulate the 
central stellar temperature. 

Unlike the accretion history, which is a macroscopic 
property that almost certainly varies from star to star, 
the initial radius is at least partly fixed by microphysics. 
The initial model should correspond to the "seed" pro- 
tostar which forms as a result of second collapse in- 
duced by collisional dissociation of hydro gen molecules 
in a t hermally-supported first core (e.g. LarsonI 119691 : 
' ^ ' 119981 : MIOO). 



iWinkle r & Newman 1980; Masunaga et al 



The entropy content of the resulting seed protostar, and 



thus its initial radius, is therefore at least partly set by 
the properties of the hydrogen molecule, in which case we 
would not expect large star-to-star variations. Nonethe- 
less, models of second collapse have not extensively ex- 
plored the influence of factors like rotation or magnetic 
fields. Thus, to be conservative we consider a factor of 
~ 100 variation in possible initial radii, extending to val- 
ues both larger and smaller than the radius of 4 Rq com- 
puted by MIOO. The smallest radii we consider are much 
smaller than have been produced in any calculation of 
second collapse. We list the full set of models we exam- 
ine in Table [TId. 

Figure[5]presents the evolution of the stellar radius un- 
til mass accretion ceases in each case. We see that for 
all stellar masses these cases vary much more than those 
in Figure [5] With the smaller initial radius, the stel- 
lar radius is also smaller after the stellar mass increases 
by mass accretion. In case mC5-C-Ri0.3, for example, 
the temperature in the stellar interior is initially higher 
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Fig. 5. — Same as Fig. EJbut for cases with diflferent initial 
radii: Rt,o = 8.0 Rq (green dotted, case mC5-C-Ri8), /?*,o = 
3.7 Rq (green dashed, case mC5-C-Ri3.7), i?*,o = 1-5 ^0 (red 
solid, case mC5-C), 1.0 Rq (red dotted, case mC5-C-Ril), 0.65 Rq 
(red dashed, case mC5-C-Ri0.65), 0.3 Rq (blue solid, case mC5- 
C-RiO.3), 0.25 Rq (blue dotted, case mC5-C-Ri0.25), and 0.2 Rq 
(blue dashed, case mC5-C-Ri0.2). The values of the initial radii 
are labeled in the panel. The cold mass accretion at the constant 
rate M = 10~^ Mq yT~^ is adopted for all the cases. 
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Fig. 6. — The evolution tracks with varying the initial radius 
with the cold mass accretion. The constant accretion rate of 
10~^ Mq yr~^ is adopted for these all cases. The tracks un- 
til the protostellar mass reaches 0.9 Mq by mass accretion are 
presented. The different lines shows the evolution with differ- 
ent initial radii: i?*,o = 8.0 -Rq (green dotted, case mC5-C-Ri8), 
R*fi = 3.7 -Rq (green dashed, case mC5-C-Ri3.7), -R*,o = 1-5 Rq 
(red solid, case mC5-C), 1.0 -Rq (red dotted, case mC5-C-Ril), 
0.65 -Rq (red dashed, case mC5-C-Ri0.65), 0.3 -Rq (blue solid, 
case mC5-C-Ri0.3), 0.25 Rq (blue dotted, case mC5-C-Ri0.25), 
and 0.2 -Rq (blue dashed, case mC5-C-Ri0.2). The initial radius 
in each model is labeled in the figure. The input parameters for 
these models are summarized in Table [l}). 



compared to case mC5-C. As a result, deuterium burning 
begins earlier at M,» ~ 0.03 Mq, and the stellar radius is 
always smaller than 1 Rq. 

Figure [7] shows the evolutionary tracks in the HRD 
until the stellar mass reaches 0.9 Mq. We see that tracks 
with small initial radii occupy the lower part of the HRD, 
reflecting the variation shown in Figure[5l Figure[7]shows 
the evolution of the stellar positions in the HRD after 
mass accretion ceases for these cases, computed for stars 
that stop accreting at final masses of 0.05 Mq, 0.1 Mq, 
0.3 Mq, 0.5 Mq, 0.7 Mq, and 0.9 Mq. In comparison 
to Figures [3] and m we see a much larger spread between 
the isochrones than that produced by different accretion 
histories. The stellar positions significantly differ even 
for the same mass and age. We see that the observational 
data points near the 10 Myr non-accreting isochrone are 
covered even in the snapshots for t < 1 Myr. 

However, with the exception of models mC5-C-Ri3.7 
and mC5-C-Ri8, which start with large initial radii and 
entropies, the distribution of the calculated PMS stars 
is never consistent with that of the observational data 
points. Even at the earliest time snapshot, in these mod- 
els most or all the stars with Tcs > 3500 K lie below 
the 10 Myr isochrone, where there are no observed stars. 
This problem is particularly serious for the cases repre- 
sented with blue symbols (cases mC5-C-Ri0.3, mC5-C- 
RiO.25, and mC5-C-Ri0.2). 

In order to render these cold accretion models with 
small radii consistent with observations, one would have 
to posit that only stars whose final masses are below 
0.5 Mq have second cores with radii much smaller than 
the values predicted by MIOO and similar calculations. In 
effect, the ~ 0.01 Mq second core would need to know 
in advance what the final properties of the star would 
be, and the properties of the second core would have to 
somehow correlate with the final mass. Given the limited 
range of second collapse models that have been explored 
in the literature, we cannot rule out the possibility that 
both very small second core radii and a systematic corre- 
lation between second core radii and final stellar masses 
exist in nature. However, we are also unaware of any 
observational or theoretical evidence in favor of either of 
these propositions. 

If, on the other hand, we restrict our attention to 
models with initial radii such that the stars are at least 
marginally consistent with the data at all effective tem- 
peratures (models mC5-CRi-8 to mC5-CRi0.65, shown 
in green and red in Figure [71 we see that the spread in 
HRD location at Tes < 3500 K is very small, comparable 
to the spread seen in Fiugres |3] and 01 and much smaller 
than the spread in observed stellar positions. 

3.3. Thermal Efficiency Variation 

Next, we consider the effects of different thermal effi- 
ciencies of mass accretion on protostellar evolution. In 
addition to the cold accretion models we have previously 
considered we now add models that experience hot ac- 
cretion for varying lengths of time during their evolu- 
tion. We summarize these models in Table [TJ;. Model 
mC5-H uses pure hot accretion, and models mC5-HCa; 
switch from hot to cold accretion once their masses ex- 
ceed X Mq. 

For these new hot models we adopt 0.01 Mq radia- 
tive stars for the initial state, as in the previous section. 
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Fig. 7. — Same as Fig. |3]but for the cases presented in Fig. [6]and for tlie PMS stars of various masses. The symbols in each panel mark 
the positions of 0.9 Mq, 0.7 Mq, 0.5 Mq, 0.3 Mq, 0.1 Mq, and 0.05 Mq stars from the left for cases mC5-C-Ri3.7 (green filled circles), 
mC5-C (red asterisks), mC5-C-Ril (red open circles), mC5-C-Ri0.65 (red filled circles), mC5-C-Ri0.3 (blue asterisks), mC5-C-Ri0.25 (blue 
open circles), and mC5-C-Ri0.2 (blue filled circles). We omit the evolutionary tracks for clarity here. We also omit model mC5-C-Ri8 
because it is nearly identical to mC5-C-Ri3.7. 



We use initial radii of 3.7 R q-, , which is roiighly c onsistent 
with the value of 3.45 Rq in lStahler et al.l (|1980[ ) and the 
value of 4.0 Rq computed by MIOO, and somewhat larger 
than the value of 2 adopted in .Stahler (,1988) F| How- 
ever, the initial radius doe s not matter in the hot case 
as it does in the cold one (.Stahlerlll988t iHartmann et alj 
|1997[) . In the cold case, the entropy of gas accreting onto 
the star matches the entropy of the stellar atmosphere. 
Since the structure of the stellar atmosphere depends on 
the initial model, so does the subsequent evolution. In 



2 The difference between Stabl er et al.l i flQSOl) and|Stahlcr (1988) 
arises because [Stabler (1988) use fully convective initial models, 
although the interior of a seed protostar would be radiative prior 
to deuterium burning (Stabler et al. 1980). In our calculations, we 
follow the early evolution, where a convective layer occurs after the 
ignition of deuterium burning. 



contrast, for hot models there is a self-regulation mecha- 
nism that removes the dependence on the initial radius. 
If the initial radius is too small, accreting gas releases a 
large amount of gravitational energy before reaching the 
stellar surface. Since a fraction of this energy is trapped 
in the accreting gas in the hot case, accreting materials 
settle onto the star with high entropy, which increases 
the stellar radius. On the other hand, if the initial ra- 
dius is too large, the opposite effect operates. Accreting 
material has less entropy because it converts less of its 
gravitational energy to kinetic energy, and this serves to 
decrease the stellar radius. The stellar radius is regulated 
as a result. 

Figure[8]shows the evolution of the stellar radius for the 
runs with varying thermal efficiency. In the cases with 
the early hot accretion, the stellar radius at M* = M*_hc 
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Fig. 8. — Same as Fig. |2] but for cases with different thermal 
efficiencies of mass accretion. Only hot and cold mass accretion is 
adopted with the same initial radius 3.7 Rq in cases mC5-H (pur- 
ple) and mC5-C-Ri3.7 (green). In the cases with the magenta lines, 
the hot mass accretion is switched to the cold accretion when the 
stellar mass reaches M,^hc = 0-3 Mq (magenta dotted, case mC5- 
HC0.3), 0.1 Mq (magenta soUd, case mC5-HC0.1), and 0.03 Mq 
(magenta dashed, case mC5-HC0.03). The values of M* jjc a^re 
labeled in the figure. The input parameters for these models are 
summarized in Table [l};. The red line represents the evolution in 
case mC5-C for a comparison. 
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Fig. 9. — Same as Fig. |6]but for varying the thermal efficiency 
of mass accretion. Different lines indicate the different thermal ef- 
ficiencies: only hot mass accretion (purple, case mC5-H), cold ac- 
cretion followed by early hot accretion for M* < M.^hc = 0-3 Mq 
(magenta dotted, case mC5-HC0.3), 0.1 Mq (magenta solid, case 
mC5-HC0.1), 0.03 Mq (magenta dashed, case mC5-HC0.03), and 
only cold mass accretion with the same initial model as in the above 
cases (green, case mC5-C-Ri3.7) 



is ~ 3 Rq, significantly larger than in any of the purely 
cold cases considered in Section [3^ In case mC5-HC0.I, 
for example, the stellar radius slightly decreases after the 
boundary condition is switched to cold mass accretion 
at = 0.1 Mq. However, deuterium burning begins 
soon after and the star expands for M^, > 0.3 Mq. The 
evolution at M, > M^^c is close to that in the purely 
hot case mC5-H. Even with our lowest switching mass 
M*,HC = 0.03 Mq (case MC5-HC0.03), the stellar radius 



always exceeds 1 Rq. 

Figure [9] is the same as Figure |6] but for the cases with 
varying efficiency. We see that at the end of accretion 
all the models with any hot accretion, even mC5-HC0.03 
where we switch to cold accretion at 0.03 Mq, lie at or 
above the 1 Myr non-accreting isochrone, and above the 
level of the observed data. The cold models are sys- 
tematically lower. Significantly, only the hot models are 
consistent with the stars at the highest observed L and 
Tcs- We therefore conclude that either these stars must 
experience some hot accretion, or that their initial radii 
must be very large, thereby achieving the same effect. 

Figure [TU] shows the subsequent time evolution of the 
stellar positions for these cases as in Figure [T) We omit 
cases mC5-HC0.1 and mC5-HC0.3 here, because they 
are almost the same as case mC5-H. The snapshot at 
t = 0.3 Myr shows that, in case mC5-H, the stars are 
above the 1 Myr isochrone when mass accretion ceases. 
The PMS stars descend in the HRD after this, and their 
positions are nearly consistent with the isochrones at 
t > I Myr. The differences in stellar positions between 
cases mC5-H and mC5-HC0.03 are small even in the early 
snapshots when t < 1 Myr, and at T^g < 3500 K the 
differences are particularly small. Thus, even a short du- 
ration of hot mass accretion produces a similar result to 
the purely hot case. BCG09 also find very bimodal out- 
comes between their hot cases {a > 0.2) and cold cases 
(a = 0). This suggests that it is difficult to freely pop- 
ulate the PMS stars in the HRD only by varying the 
accretion thermal efficiency. 

4. THE RELIABILITY OF NON-ACCRETING ISOCHRONE 
ESTIMATES OF AGES AND AGE SPREADS 

Having explored the parameter space thoroughly, we 
are now in a position to make some general statements 
about the reliability of theoretical isochrones as tools for 
estimating ages and age spreads. In this analysis, it is 
helpful to separate the cases of stars with T^g > 3500 
K from those with T^g < 3500, because the results are 
quite different in the two cases. 

At Tcff > 3500 K, stars with high and low thermal ef- 
ficiencies, or with different initial radii, can end up quite 
far apart in the HRD even at equal ages (Figures [7] and 
fTU]) . We therefore conclude that, in this regime, stellar 
age and mass alone cannot uniquely determine stellar po- 
sitions in the HRD, and stellar age estimates based on 
HRD positions cannot in general be considered reliable. 
However, we note that errors in this regime only occur in 
one direction: young stars can appear old, but old stars 
never appear young. In no case do we find stars above 
the 1 Myr isochrone whose ages are actually > 1 Myr. 
Thus very young age estimates are reliable in this effec- 
tive temperature range, even if old age estimates are not. 
With regard to age spreads, we note that the snapshots 
at t < 1 Myr show that the models in which we vary the 
thermal efficiency and the initial radius cover the entire 
observed spread of PMS stars with Toff > 3500 K. Thus, 
the entire observed luminosity spread in this temperature 
range could be explained if a coeval population were to 
consist of some stars that underwent hot accretion and 
others that underwent cold accretion. Estimates of age 
spreads in this regime are therefore unreliable. 

The situation is quite different for T^g < 3500 K. In 
this regime, consulting Figures [3j lU [71 and [10] shows 
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Fig. 10. — Same as Fig. [71 but for the cases presented in Fig. (9] The symbols in each panel mark the positions of 0.9 Mq, 0.7 AIq, 
0.5 Mq, 0.3 Mq, 0.1 Mq, and 0.05 Mq stars from the left for cases mC5-H (purple filled circles), mC5-HC0.03 (magenta filled circles), 
mC5-C-Ri3.7 (green filled circles), and mC5-C (red asterisks). 



that the only models capable of producing false old ages 
are those that have pure cold accretion starting from 
very small initial radii. For example, even at a true 
age of 3 Myr, only model mC5-C-Ri0.2 places stars with 
Tcff = 3000 K below the 10 Myr isochrone. However, 
these models can be considered acceptable and consis- 
tent with the observational data only if two uncertain 
propositions hold. First, a significant number of stars 
would have to accrete from seeds whose radii are an or- 
der of magnitude smaller than any that have been pro- 
duced in theoretical models thus far. Second, in order to 
avoid drastically overpopulating the region below the 10 
Myr isochrone at Teff > 3500 K, such small initial radii 



would have to be realized only for stars that have small 
final masses at the end of accretion. There would have 
to be some mechanism to generate a correlation between 
the radii of ~ 0.01 Mq second cores and the masses of 
the ~ 1 Mq stars to which they eventually grow. Given 
the limitations of the theoretical models, neither of these 
propositions can be ruled out, but there is currently no 
evidence in favor of them either. 

If we rule out the cold accretion, small initial radius 
models based on these problems, we find that the ob- 
served luminosity spread for Toff < 3500 K is consider- 
ably larger than the spread among the remaining mod- 
els. All these models lie fairly close to the non-accreting 
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isochrones appropriate to their true ages for times t > 3 
Myr, and at times t > 1 Myr at T^s - 3000 K. We 
therefore conclude that, unless ^ 0.01 Mq seed proto- 
stars have some very specific unexpected properties, age 
estimates of 1 — 3 Myr or more based on non-accreting 
isochrones are reliable for stars with effective tempera- 
tures below ~ 3500 K, at least insofar as the observations 
themselves are reliable. Neither variation in thermal effi- 
ciency, accretion history, or initial radius can explain the 
observed spread in the HRD. 

5. SUMMARY AND DISCUSSION 

In this paper, we have examined a variety of low- 
mass protostellar evolutionary tracks with varying ac- 
cretion histories, initial models, and thermal efficiencies 
of mass accretion. We have also compared the resul- 
tant spread of PMS stars in the HRD t o that observed 
in nearby low-mass sta r forming regions (IPeterson et al.l 
l2008t IGatti et al.ll2Q0l 2008; M uzeroUe et al.ll20Q5D . 

We first calculate protostellar evolution models with 
varying accretion histories but fixed initial stellar mod- 
els and boundary conditions (Sec. 13. 1() . Our results show 
that if mass accretion is thermally inefficient, variation 
in the accretion history hardly infiuences protostellar 
evolution. Although isochrones f or non-accreting proto- 
stars, s uch as those ca l culated bvlD'Antona fc Mazzitellil 
(fTOQl : IBaraffe et all (IT9981) : ISiess et al.l (120001) . do not 
necessarily provide us with correct stellar ages, models 
with differing accretion histories nevertheless form a tight 
sequence in the HRD. Thus variable accretion histories 
alone cannot explain the observed spread of PMS stars 
in the HR. Moreover, the errors in absolute ages arise be- 
cause non-accreting isochrones are not good descriptions 
of stars growing with thermally inefficient mass accre- 
tion. They are not a result of variable accretion histo- 
ries. However, we note that this does suggest that accret- 
ing isochrones may resolve the problem of systematically 
larger inferred ages for high-mass clu ster members (e.g., 
IHillenbrandl[200a ICovev eTalllMoh . 

Second, we examine protostellar evolution with differ- 
ent initial models and thermal efficieiicies of mass accre- 
tion, using a constant accretion rate M = 10~^ Mq yr~^ 
(Sec.OandESI). We find that the spread of PMS stars 
in the HRD that results from varying the initial radius 
(or entropy) or the thermal efficiency is much larger than 
the spread that arises from different accretion histories. 
We find that a coeval population of stars with signifi- 
cant star-to-star variation in thermal efficiency or initial 
radius could conceivably occupy the entire observed lu- 
minosity range for protostars with effective temperatures 
> 3500 K. Thus ages and age spreads observed in this 
temperature range may be unreliable. At lower effec- 
tive temperature, however, the situation is very differ- 
ent. The only models we found in our parameter space 
survey that are capable of producing false old ages at 
low Teff are those with purely cold accretion starting 
from very small initial radii. However, we point out that 
these models require that 0.01 Af© second protostellar 
cores have radii < 0.3 i?©, more than an order of mag- 
nitude smaller than any formed in simulations to date 
(e.g., MIOO), and that they can be rendered consistent 
with observations only if such small initial radii are re- 
alized only for stars that end up growing to small final 
masses. Neither possibility can be definitively ruled out. 



but neither is supported by any current observations or 
theory either. If we exclude very small initial radii on 
these grounds, we find that all the remaining models 
indicate that ages and age spreads inferred from non- 
accreting isochrones are reliable for cool stars, at least 
to the extent that the observationally-det ermined lumi- 
nositie s and temperatures are reliable (e.g.. lDa Rio et all 

i2oTobe . 

In varying the thermal efficiency, we find that models 
with only a small amount of hot accretion, e.g., while 
M* < 0.03 Mq, nonetheless show similar evolution to 
models with entirely hot accretion. Thus, in order to 
explain the HRD spread at low-masses with cold ac- 
cretion models, the models must begin with small radii 
and be thermally inefficient for nearly their entire evolu- 
tion. Observational constraints on the thermal efficiency 
of accretion are somewhat tenuous. During the earliest 
stages during which accretion rates of 10~^ — 10"'* M© 
yr~* occur, protostars are too deeply embedded with 
too much radiation reprocessing to measure accretion 
signatures directly. Observations of T Tauri stars, for 
which accretion rates have declined to < 10~^ M© yr~^, 
suggest that for these rates the accretion column cov- 
ers only ~1-10% of the stellar surface. This supports 
a cold accretion senario later in the accretion history. 
However, observations also find t hat the covering frac- 
tion increas es with accretion rate ([Gullbring et al.ll2000l : 
lArdila fc^ asri 2000), suggesting that accretion may be 
more thermally inefficient at early times. Higher mass 
stars, which experience higher initial accretion rates, 
seem unlikely to avoid hot accretion, while it is more 
observationally probable that very low-mass stars expe- 
rience purely cold accretion. 

Our conclusions are different from those of BCG09, 
who stressed the significance of episodic mass accretion 
for explaining the observed spread of PMS stars in the 
HRD. However, our numerical results are actually con- 
sistent with theirs. First, BCG09 calculated protostellar 
evolution using simple, non-episodic accretion histories 
(their Fig. 1) and found that, with thermally efficient 
accretion (a > 0.2 in their notation), isochrones for non- 
accreting protostars give the correct ages at i > 1 Myr. 
Disagreement with the isochrones arises only for ther- 
mally inefficient accretion flows (a = 0). This is consis- 
tent with the results shown in our Figure [T) Next, they 
calculated the evolution with more vigorous, episodic 
mass accretion histories (their Fig. 2). However, the 
spread they obtain is no broader than that shown in their 
Figure 1, indicating that episodic accretion does not in- 
crease the HRD spread beyond what they had already in- 
troduced by using varying initial conditions and thermal 
efficiencies. Indeed, they state that the time-dependence 
of the accretion rates is not essential for their results, 
and we confirm this finding^ Our results suggest that 
BCG09 were able to obtain small luminosities for stars 
with Tcs < 3500 K, and thus claim to reproduce the 
observed HRD using a coeval population, because they 
used cold accretion starting from initial conditions with 
entropies far lower than what current theoretical mod- 

^ IBaraffe fc Chabrieii II2010I ) argue that episodic mass accretion 
also significantly influences the lithium depletion of low-mass stars. 
However, our present work suggests that variation of the initial 
models will be more significant than variation of the accretion his- 
tories for the problem of lithium depletion as well. 
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els predict. They also did not continue runs with these 
initial conditions up to higher masses and effective tem- 
peratures|3 Our results suggest that, had they done so, 
the resulting stars would have fallen well below the locus 
of observed stars in the HRD, as do our comparable low 
initial entropy models. 

Finally, we stress that we do not reject the episodic 
mass accretion as a possible solution for the "lu minos- 
ity problem" of young ei i ibedd e d sources ( Dunhair Tet al.l 
l2nTnt IMcKee fc OfaMl IMfil : lOffner fc McKeel 1201 Ih . 
Episodic accretion may well occur. It is simply not ca- 
pable of explaining the broad spread of optically-visible 
PMS stars in the HRD. 
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* For example, the only cases BCG09 show where stars with 
Tcfl < 3500 K lie near the 10 Myr non-accreting isochrone after 
1 Myr are their models A-C. However, they only use the initial 
conditions and accretion rates for models A-C to produce stars up 
to 0.2 Mq. In comparison, all the models that they run to masses 
above 0.5 Mq have much higher starting entropies. 
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